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ABSTRACT 

We study the numerical convergence of hydrodynamical simulations of self-gravitating 
accretion discs, in which a simple cooling law is balanced by shock heating. It is well- 
known that there exists a critical cooling time scale for which shock heating can no 
longer compensate for the energy losses, at which point the disc fragments. The nu- 
merical convergence of previous results of this critical cooling time scale was questioned 
O^l recently using Smoothed Particle Hydrodynamics (SPH). We employ a two-dimensional 

grid-based code to study this problem, and find that for smooth initial conditions, frag- 
mentation is possible for slower cooling as the resolution is increased, in agreement with 
p j recent SPH results. We show that this non-convergence is at least partly due to the 

PjJ creation of a special location in the disc, the boundary between the turbulent and the 

laminar region, when cooling towards a gravito-turbulent state. Converged results appear 
to be obtained in setups where no such sharp edges appear, and we then find a critical 
cooling time scale of ~ 4S1~ 1 , where fl is the local angular velocity. 

Key words: planets and satellites: formation -planetary systems: protoplanctary discs 
- accretion discs - hydrodynamics - instabilities 
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1 INTRODUCTION 

In the absence of any heating mechanism, a self-gravitating 
accretion disc will cool down until it becomes gravitation- 
ally unstable, which happens when the stability parameter Q 
( |Toomre|1964] > 

Here, c s is the sound speed, k denotes the epicyclic frequency, 
and E is the surface density of the disc. When Q approaches 
unity and gravitational instabilities kick in, spiral shocks are 



generated that heat up the disc ( Goldreich fc Lynden-Bell 
|1965 Durisen et al.|200~ \ . An equilibrium can then be set up 
in which shock heating exactly balances energy losses through 
cooling ( Paczynski|1978 1 . This is called the gravito-turbulent 
state of the disc. If the cooling time scale is parametrised as 



^cool — p\\ 



(2) 



with f3 constant and Q the local angular velocity of the disc, 
a local balance of heating and cooling dictates that the spi- 



ral shocks generate an effective a parameter (Pringle 1981 
Gammie||2001[ ) 
1 



4 

a = - 



7(7-1)/?' 



(3) 



where 7 is the adiabatic exponent. A lot of work has been put 
into determining whether a local equilibrium really is estab- 
lished and whether an a-description is valid at all for grav- 



2011 ), which is not necessarily the case (Balbus & Papaloizou 



1999). 



itoturbulent discs (Lodato & Rice 2004 2005 Forgan et al 



For strong enough cooling, the disc is unable to provide 
enough shock heating to balance the cooling, and the disc 



fragments (Gammie 20011, possibly forming gas giant plan- 
ets. The critical value of P was found to be f3 c « 3 i n|Gammie| 



(20011 for 2D simulations with 7 = 2. Rice et al. (20051 in 



terpreted the fragmentation criterion as a maximum stress 
the disc can sustain. They found a max ~ 0.06, which leads 
to a value of (5 C that depends on 7 through equation jjjl. It 
was found by Clarke et al. ( 2007 I that the value of a ma x can 
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increase by a factor of ~ 2 when /3 is decreased slowly. 

The determination of the critical value of /3 is important 
in models of giant planet formation that rely on disc insta- 
bility (e.g. |Boss|[i997[ |Boley et al.||2010[ ). Since according to 
the results above, the cooling time scale should be of the or- 
der of the orbital time scale, planet formation through disc 
instability is more likely to occur in the outer parts of real 
protoplanetary discs ( |Rafikov|2005"} |Boley||2009[ ) . 

More recently, the numerical convergence of the results 
on /3 C was questioned by Meru & Bate ( 2011a|b ' I. There, it 
was shown through SPH simulations that the critical value 
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of P had not converged with respect to resolution. At higher 
and higher resolution, discs with higher and higher j3 would 
fragment, with no sign of convergence. It was suggested by 
Lodato & Clarke ( 20111) that this behaviour might be due 



to the SPH artificial viscosity or the artificial smoothing of 
density enhancements. 

In this Letter, we present results on the convergence of 
the critical value of f3 using the grid-based hydrodynamics 
code FARGO (Fast Advection in Rotating Gaseous Objects, 
Masset|2000a|b see also section|2|. First of all, with a smooth 
initial setup as in many previous simulations of gravitationally 
unstable discs, we find similar non-convergence with FARGO, 
indicating that the problem is not specific to SPH. We then 
show that numerical convergence can be reached by avoiding 
a sharp boundary between the gravito-turbulent inner part of 
the disc and the still laminar outer part, which arises when 
using smooth initial conditions during the initial stage of cool- 
ing towards Q ~ 1. If smooth initial conditions are not used, 



we find a value of /3 C close to that found by Gammie (2001 1. 



2 NUMERICAL METHOD 



We use the grid-based hydrodynamics code FARGO (Masset 
2000a|b|), for which a self-gravity solver was presented in 
Baruteau & Masset (20081. We have implemented the sim 



pie cooling law 



de 
dt 



t 

^cool 



(4) 



where e is the internal energy and t coo i is given by equation 
(|2[). A standard prescription for artificial viscosity (Stone & 



Norman |1992| is used to handle shocks. Following Lodato 
fc Rice| p004| , we choose our disc to lie between R = 0.25 
and R = 25. The self-gravity module of FARGO requires a 
logarithmic grid in the radial direction, and we choose our 
grid so as to give square cells everywhere (AR/R « A<f>). Our 
lowest resolution has 512 cells in the radial direction and 768 
cells in the azimuthal direction, and we increase the resolution 
by a factor of 2 and 4. We use outflow boundary conditions 
at the inner and outer grid edges. 



Following Lodato & Rice ( 2004 I , we choose our disc to 
have surface mass density E oc R~ and temperature T oc 
rp^g an g U i ar velocity is Keplerian, corrected for the 
mass and pressure of the disc. A small level (~ 0.1%) of white 
noise is put on top to break the axisymmetry and allow spiral 
waves to form. The surface density and aspect ratio H/R at 
R = 1 are chosen so as to give a total disc mass of 0.25 M*, 
and so that the minimum value of Q ~ 2 at the outer edge of 
the disc. We use 7 = 5/3 throughout. 

We measure the total shear stress in the simulations by 
calculating 



Tru: 



gag-? 

4nG' 



dz \ + {Y,8vr8v v ) . 



(5) 



where gu and g v are the gravitational accelerations, 8vr and 
8v v are velocity fluctuations, and (•} denotes an azimuthal 
average. The vertical integral is calculated by computing gn 



and g v at different values of z (see the Appendix of Baruteau 
et al.|2011[ for details). To compare with equation ([3jl, we can 
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Figure 2. Radial profile of Q (top panel), H/R (middle panel) 
and S (bottom panel) for a resolution AR/R = 0.008 after (solid 
curve), 50 (dotted curve), 100 (dashed curve) and 150 (dash-dotted 
curve) orbits at R = 1 for j3 = 7.5. 



use an a-parametrisation: 



dlogfiX Tj 



'fly 

dlogRj Ec s 2 ' 



(6) 



We usually find that the disc settles into a state of constant 
a, in which case a single value of a can be assigned to the 
disc. 



3 RESULTS 

At first, the disc is almost axisymmetric with no source of 
heating. It will therefore start to cool at constant surface 
density towards Q = 1. Since j3 is constant, the cooling time 
scale icooi is shortest in the inner parts of the disc, and this 
is where Q = 1 is reached first. Therefore, gravito-turbulence 
starts from the inside and moves outward as more and more of 
the disc cools towards Q = 1. At this point, the initial tem- 
perature distribution has been washed out completely, and 
the temperature is now given by the requirement of Q w 1 
with the initial surface density. It is easy to see that this gives 
H/R oc R for our initial surface density profile. The number 
of grid cells per scale height is therefore an increasing func- 
tion of R, which means that the outer parts of the disc are 
better resolved. The disc then starts to evolve viscously to- 
wards a steady-state accreting solution with constant values 



of Q and a, which has E oc R 3 ^ 2 and constant temperature 
(see |Pringlef[98l) equation (2.10)). 

We choose /3 = 7.5 and evolve the discs for 2000 or- 
bits at R = 1. According to Rice et al. ( 2005[ ) and |Gammie 
( |2001| ), these discs should not fragment for 7 = 5/3, though 



previous results may be affected by resolution ( Meru & Bate 



2011a I. We increase the resolution to see if convergence can 



be reached. The resulting surface density distributions are 
shown in Fig. [I] In our lowest resolution case (left panels), 
we resolve H outside R « 3, but we find no evidence of frag- 
mentation (defined such that the fragment is more than 2 
orders of magnitude denser than the surroundings). In the fi- 
nal gravito-turbulent state (bottom left panel), we measure a 
total stress that compares very well with the predicted value 
of equation |3| . 

When increasing the resolution by a factor of 2 (middle 
panels of Fig. |T|) , we do find fragmentation. After 150 orbits at 
R — 1, clumps start to form at the outer edge of the turbulent 
region (top middle panel in Fig.[IJ). The first surviving clumps 
form at R ~ 15, even though the necessary length scale (H) is 
resolved much further in. When increasing the resolution even 
further (right panels of Fig. [I]), the disc fragments further in 
(R ~ 6). After 300 orbits at R = 1, several of the fragments 
have merged or migrated off the computational domain (see 



Baruteau et al. 12011 1, leaving five distinct clumps. 



We have observed fragmentation for values of j3 as high 
as 15, at a resolution that is 8 times the lowest resolution 
presented in this Letter. Interestingly, for /3 = 15, clumps 
form in the outer disc, migrate inward and leave the com- 
putational domain. After all clumps have left the disc, no 
further fragmentation is observed. This already suggests that 
the fragmentation for high ft is transient and related to the 
initial state of the disc. 

We can obtain further insight into why these discs frag- 
ment by looking at radial profiles of Q, H/R and E. Figure 
[2] shows the results for the low resolution simulation with 
/3 = 7.5. From the top panel, we see that the disc cools down 
to a state of constant Q « 1.5. Note that inside R = 3, the 
scale height the disc tries to cool towards can not be resolved 
at this resolution, so we should only consider R > 3. As men- 
tioned before, the disc becomes turbulent from the inside out. 
The outer edge of the turbulent region can be clearly iden- 
tified as a location where Q < 1 in the top panel, and as 
bumps in the temperature and surface density profiles in the 
bottom two panels. Such a global feature makes the set-up 
no longer scale- free. A region where Q < 1 is of course unsta- 
ble to gravitational instabilities. Moreover, sharp transitions 
in temperature and surface density are notoriously unstable 
(e.g. |Papaloizou fc Pringle||1985| p87] |Lovelace et al.||1999 



Lin & Papaloizou 2011 1. It is conceivable that the formation of 



clumps at the outer edge of the turbulent region, as observed 
for example in the top middle panel of Fig. [TJ is related to 
this sharp transition between the turbulent and laminar parts 
of the disc. 

To test to what extent the smooth initial conditions, and 
the related formation of this global feature, play a part in 
disc fragmentation, we employ the same strategy as |Clarke| 
et al. (20071. We start off with a disc that has /3 — 30, and let 



gravito-turbulence fully develop over 1000 orbits at R = 1. 
For this value of /3, we have not observed fragmentation for 
any of the resolutions. Over the next 500 orbits, we then lin- 
early decrease /3 to a value of 15, 7.5, 5, 4 or 2.5. This way, we 
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Figure 3. Surface density after 2000 orbits at iZ = 1, where in the first 1000 orbits f3 = 30 was used, after which was decreased to 7.5 in 
the next 500 orbits. From left to right, the resolution is AR/R = 0.008, 0.004 and 0.002. 
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Figure 4. Measured total shear stress for different final values of /3 
at different resolutions. The value of a was calculated by averaging 
over orbits 1900 — 2000 of each run (i.e. when the discs were held 
at the final constant /3 values), and averaged over the whole disc. 
The dotted line indicates equation |3j. 
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avoid any initial transients affecting the possible fragmenta- 
tion of the disc, and restore the scale- free nature of the set-up. 
We then evolve the discs with /3 held fixed at the desired val- 
ues for another 2000 orbits (or until they fragment). 

In Fig. [3J we show the resulting surface density struc- 
tures after 2000 orbits for the same three resolutions and the 
same final /3 = 7.5 as in Fig. [1] This time, we do not observe 
fragmentation for any of the resolutions, and the surface den- 
sities appear very similar. We also find that the resulting to- 
tal stresses are equivalent in the three cases (see Fig. [4| , and 
agree very well with equation (j^3j) . We therefore conclude that 
we have reached numerical convergence in this case. A side- 
effect of letting the disc evolve for 3500 orbits is that mass is 
redistributed over the computational domain, as the disc is 
trying to reach a steady state with £ oc R~ 3 ^ 2 . 

We find that /3 can be reduced to a value of 5 without 
seeing fragmentation for all 3 resolutions. From Fig. [4] we see 
that in all cases, we find very similar total shear stresses for 
different resolutions, which in turn agree very well with equa- 
tion J3j. For p < 4, we find fragmentation for all resolutions 
considered. The maximum stress the disc can provide to bal- 
ance the cooling is a max « 0.1. This is in agreement with the 



results of Clarke et al. (20071, who found that the maximum 



Figure 5. Surface density after 1550 orbits at R = 1, of which the 
first 1000 orbits were run with /3 = 30, and in the next 500 orbits 
/3 was decreased to 2.5. The disc fragmented after 1475 orbits. 
Resolution is AR/R = 0.004. 



Rice et al. (2005), who started from smooth initial conditions. 



Note that we do not attempt to determine the exact value of 
/3 C - Around f3 = /3 C , other numerical effects may play a role 
(see |Clarke et aL|2007[ ) . 

Figure [5] shows the fragmented disc with /3 = 2.5 at 
AR/R — 0.004. Note that there are no special locations any- 
more: the disc fragments everywhere as long as H is resolved, 
which is basically the criterion of Truelove et al. (19971. We 



also do not find any borderline cases ( Meru & Bate 2011a I, for 



stress increased by approximately a factor of 2 compared to 



which the disc fragments at early times but the fragments do 
not survive. The disc either fragments, after which the formed 
clumps migrate and merge, or the disc finds itself in a steady 
gravito-turbulent state. With this setup, a single parameter, 
f3, determines whether the disc fragments or not. 
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4 DISCUSSION AND CONCLUSIONS 



ACKNOWLEDGEMENTS 



We have studied the numerical convergence of simulations of 
self-gravitating discs using the grid-based code FARGO. For 
smooth initial conditions, we find the same non-convergence 
as recently reported for SPH simulations (Meru & Bate 



|2011a[ ). We have argued that this non-convergence is related 
to the smooth initial conditions, that lead to an edge in tem- 
perature and surface density at the outer boundary of the 
gravito-turbulent region. This global feature, which becomes 
more pronounced at higher resolution, can drive instability at 
larger values of /3, leading to fragmentation at least for /3 ^ 15. 
It is not clear whether convergence can ever be reached using 
this set-up. It is beyond the scope of this Letter to characterise 
these edge effects in detail, because the edge is completely 
artificial. We comment that for a constant initial surface den- 
sity, the disc fragments for even higher values of (5 at a given 
resolution (/3 = 15 for AR/ R = 0.002, but note that, since 
H/R oc R 2 in this case, the number of grid cells per scale 
height increases rapidly towards the outer disc). 

It is possible to restore the scale-free nature of the set up 
by carefully choosing the initial conditions. We have taken the 
approach of |Clarke et al.| ( |2007[ ), and first set up a gravito- 
turbulent disc at a high value of /3 = 30, for which we do not 
find fragmentation with the resolutions we consider. We then 
linearly decrease /3 to the desired value. Using this set up, 
we find that our results are converged numerically, and the 
maximum value of a that the disc can sustain, a max ~ 0.1, is 
in good agreement with that fo und in|Clarke et al.| pb07| ). It 
is important to note that while Clarke et al. ( 20071) interpret 



their results as being partly dependent on the disc's thermal 
history, their results may well have been due to the removal 
of edge effects. We find a critical value of f3, /3 C « 4, in reason- 



able agreement with Gammie (2001 1. Note also that Gammie 
(20011 did not find fragmentation for /3 = 10 for any of the 



resolutions he considered, which is not surprising in the light 
of our findings, since he considered local models, from which 
edges are absent by definition. Note also that at the maxi- 
mum resolution considered by |Gammie| ( |200lj ), H is still only 
resolved by ~ 10 grid cells. More work is necessary to con- 
firm numerical convergence at even higher resolution, in local 
as well as global models, where global models should slowly 
reduce /3 to avoid artefacts from the initial conditions. It is 
conceivable that at higher resolution, an even slower reduc- 
tion rate for j3 is necessary to avoid these. Reducing ft from 30 
to 7.5 within one orbit at R = 1 again leads to the formation 
of an edge and fragmentation similar to the top panels in Fig. 

We stress at this point that it may be possible for a disc 
with j3 < /3 C to resist fragmentation (e.g. if other heating 
mechanisms are present). Likewise, it is interesting to note 
that if a physical edge exists (rather than the artificial numer- 
ical one created by smooth initial conditions), it may drive a 
disc with P > /3 C towards fragmentation. In this respect, it is 
not clear what the 'correct' initial conditions are, especially 
since the self-gravitating state occurs shortly after the for- 
mation of the disc. Nevertheless, in the idealised case of no 
edges and a simple cooling law, we have shown that at least 
spurious fragmentation at high values of ft can be avoided, 
and that this result hold for all resolutions considered here. 
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